// Solve the Momentum equation

double blend = 0;
double arg, num, denom;
vector cc;
forAll(spongeFriction.internalField(), celli)
{
	cc = mesh.C().internalField()[celli];
		
	arg = 2.0*double(cc.component(0)-spongeOnset)/spongeLength;
	arg = std::max(std::min(arg, 20.0), -20.0);
	//Info << "arg "<<"is "		<<arg		<<"\n";
	num = (std::exp(arg)-1);
	//Info << "num is "		<<"is "		<<num		<<"\n";
	denom = (1e-10+std::exp(arg)+1);
	//Info << "denom is "		<<"is "				<<denom		<<"\n";
	
	blend = (.5+.5*num/denom);
		
	spongeFriction.internalField()[celli] = blend*spongeFrictionCoefficient*(U.internalField()[celli]-spongeVelocity);
}

tmp<fvVectorMatrix> UEqn
(
    fvm::ddt(rho, U)
  + fvm::div(phi, U)
  - fvm::laplacian((MU+MU_eddy), U)
  //==
  //-spongeFriction
);

UEqn().relax();

volScalarField rAU(1.0/UEqn().A());

if (pimple.momentumPredictor())
{
    solve(UEqn() == -fvc::grad(p));
}
else
{
    U = rAU*(UEqn().H() - fvc::grad(p));
    U.correctBoundaryConditions();
}
